clear; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NxN_decomposition_GZ.m
% 
% Code for "Distortions in Production Networks"
% Bigio and La'O 2020
% produces figures and tables of model implications
%
% inputs IOvars.mat, IOparams.mat, IOseries.mat, IOshock.mat
% uses network.m function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('IOvars.mat');
load('IOparams.mat');
load('IOseries.mat');
load('IOshock.mat');

global NsecGZ

beg_year=1997;              % first year
end_year=2014;              % final year
N_years=18;                 % Number of Years
NsecGZ=size(phi_GZ_mat,1);  % Number of Industries
dYearsvec=beg_year+1:1:end_year;
Yearsvec=beg_year:end_year;

onevec      =ones(NsecGZ,1);


%% Set Parameters
v_vec       =v_i_vec;
alpha_vec   =alpha_i_vec;
eta_vec     =eta_i_vec;


% Define Horizontal economy
alphahor        =.9999*onevec;
Whor            =(1/NsecGZ)*ones(NsecGZ,NsecGZ);


% symmetric shocks time series: uses average phi_GZ
xsec_mlogphiGZ=ones(N_years,1);
for jj=1:N_years
    xsec_mlogphiGZ(jj)=mean(logphi_GZ_mat(:,jj));
end

logphiGZ_sym_mat=ones(NsecGZ,1)*xsec_mlogphiGZ';
phiGZ_sym_mat=exp(logphiGZ_sym_mat);
Dlogphi_GZsym=diff(xsec_mlogphiGZ);


%% Run Simulations
% Declare Variables

% Baseline with GZ shocks
log_A_GZ_DRS                =ones(1,N_years);
neglogLambda_GZ_DRS         =ones(1,N_years);
etabarGZ_DRS                =ones(1,N_years);
% Baseline with Symmetric GZ Shocks
log_A_GZsym_DRS             =ones(1,N_years);
neglogLambda_GZsym_DRS      =ones(1,N_years);
etabar_GZsym_DRS            =ones(1,N_years);
% Horizontal with symmetric GZ shocks
log_A_hor_GZsym_DRS         =ones(1,N_years);
neglogLambda_hor_GZsym_DRS  =ones(1,N_years);
etabar_hor_GZsym_DRS        =ones(1,N_years);
% Horizontal with asymmetric GZ shocks
log_A_hor_GZasym_DRS        =ones(1,N_years);
neglogLambda_hor_GZasym_DRS =ones(1,N_years);
etabar_hor_GZasym_DRS       =ones(1,N_years);


% CRS Baseline with GZ shocks
log_A_GZ_CRS                =ones(1,N_years);
neglogLambda_GZ_CRS         =ones(1,N_years);
etabarGZ_CRS                =ones(1,N_years);
% CRS Baseline with Symmetric GZ Shocks
log_A_GZsym_CRS             =ones(1,N_years);
neglogLambda_GZsym_CRS      =ones(1,N_years);
etabar_GZsym_CRS            =ones(1,N_years);
% CRS Horizontal with symmetric GZ shocks
log_A_hor_GZsym_CRS         =ones(1,N_years);
neglogLambda_hor_GZsym_CRS  =ones(1,N_years);
etabar_hor_GZsym_CRS        =ones(1,N_years);
% CRS Horizontal with asymmetric GZ shocks
log_A_hor_GZasym_CRS        =ones(1,N_years);
neglogLambda_hor_GZasym_CRS =ones(1,N_years);
etabar_hor_GZasym_CRS       =ones(1,N_years);


Wii         =W_2007;        %Set input-output matrix to 2007 tables
for ii=1:N_years
       
    % DRS with GZ Shocks
    [log_A_GZ_DRS(1,ii),neglogLambda_GZ_DRS(1,ii),etabarGZ_DRS(1,ii)]=network(phi_GZ_mat(:,ii),logphi_GZ_mat(:,ii),Wii,v_vec,alpha_vec,eta_vec);
    % DRS with Symmetric GZ Shocks
    [log_A_GZsym_DRS(1,ii),neglogLambda_GZsym_DRS(1,ii),etabar_GZsym_DRS(1,ii)]=network(phiGZ_sym_mat(:,ii),logphiGZ_sym_mat(:,ii),Wii,v_vec,alpha_vec,eta_vec);
    % DRS Horizontal with Symmetric GZ Shocks
    [log_A_hor_GZsym_DRS(1,ii),neglogLambda_hor_GZsym_DRS(1,ii),etabar_hor_GZsym_DRS(1,ii)]=network(phiGZ_sym_mat(:,ii),logphiGZ_sym_mat(:,ii),Whor,v_vec,alphahor,eta_vec);
    % DRS Horizontal with GZ shocks
    [log_A_hor_GZasym_DRS(1,ii),neglogLambda_hor_GZasym_DRS(1,ii),etabar_hor_GZasym_DRS(1,ii)]=network(phi_GZ_mat(:,ii),logphi_GZ_mat(:,ii),Whor,v_vec,alphahor,eta_vec);
    
    
    % CRS with GZ Shocks
    [log_A_GZ_CRS(1,ii),neglogLambda_GZ_CRS(1,ii),etabarGZ_CRS(1,ii)]=network(phi_GZ_mat(:,ii),logphi_GZ_mat(:,ii),Wii,v_vec,alpha_vec,eta_CRS_vec);
    % CRS with Symmetric GZ Shocks
    [log_A_GZsym_CRS(1,ii),neglogLambda_GZsym_CRS(1,ii),etabar_GZsym_CRS(1,ii)]=network(phiGZ_sym_mat(:,ii),logphiGZ_sym_mat(:,ii),Wii,v_vec,alpha_vec,eta_CRS_vec);
    % CRS Horizontal with Symmetric GZ Shocks
    [log_A_hor_GZsym_CRS(1,ii),neglogLambda_hor_GZsym_CRS(1,ii),etabar_hor_GZsym_CRS(1,ii)]=network(phiGZ_sym_mat(:,ii),logphiGZ_sym_mat(:,ii),Whor,v_vec,alphahor,eta_CRS_vec);
    % CRS Horizontal with GZ shocks
    [log_A_hor_GZasym_CRS(1,ii),neglogLambda_hor_GZasym_CRS(1,ii),etabar_hor_GZasym_CRS(1,ii)]=network(phi_GZ_mat(:,ii),logphi_GZ_mat(:,ii),Whor,v_vec,alphahor,eta_CRS_vec);
    
end


etabarDRS=mean(etabarGZ_DRS);
etabarCRS=mean(etabarGZ_CRS);


%% Labor Wedge Calculations

% DRS
logLWedgeGZ           = -neglogLambda_GZ_DRS;
logLWedge_GZsym       = -neglogLambda_GZsym_DRS;
logLWedge_hor_GZsym   = -neglogLambda_hor_GZsym_DRS;
logLWedge_hor_GZasym  = -neglogLambda_hor_GZasym_DRS;

% CRS
logLWedgeGZ_CRS           = -neglogLambda_GZ_CRS;
logLWedge_GZsym_CRS       = -neglogLambda_GZsym_CRS;
logLWedge_hor_GZsym_CRS   = -neglogLambda_hor_GZsym_CRS;
logLWedge_hor_GZasym_CRS  = -neglogLambda_hor_GZasym_CRS;

%% Take differences

% TFP: DRS Baseline with GZ
DlogTFP_GZ_DRS              =diff(log_A_GZ_DRS);
DlogTFP_GZsym_DRS           =diff(log_A_GZsym_DRS);
DlogTFP_hor_GZsym_DRS       =diff(log_A_hor_GZsym_DRS);
DlogTFP_hor_GZasym_DRS      =diff(log_A_hor_GZasym_DRS);

% TFP: CRS Baseline with GZ
DlogTFP_GZ_CRS              =diff(log_A_GZ_CRS);
DlogTFP_GZsym_CRS           =diff(log_A_GZsym_CRS);
DlogTFP_hor_GZsym_CRS       =diff(log_A_hor_GZsym_CRS);
DlogTFP_hor_GZasym_CRS      =diff(log_A_hor_GZasym_CRS);

% Labor Wedge: DRS Baseline with GZ Shocks
DlogLWedge_GZ_DRS           = diff(logLWedgeGZ);
DlogLWedge_GZsym_DRS        = diff(logLWedge_GZsym);
DlogLWedge_hor_GZsym_DRS    = diff(logLWedge_hor_GZsym);
DlogLWedge_hor_GZasym_DRS   = diff(logLWedge_hor_GZasym);

% Labor Wedge: CRS Baseline with GZ Shocks
DlogLWedge_GZ_CRS           = diff(logLWedgeGZ_CRS);
DlogLWedge_GZsym_CRS        = diff(logLWedge_GZsym_CRS);
DlogLWedge_hor_GZsym_CRS    = diff(logLWedge_hor_GZsym_CRS);
DlogLWedge_hor_GZasym_CRS   = diff(logLWedge_hor_GZasym_CRS);


%% Figures for Decomposition of Network Effect
% Figures in main text for network multiplier


figure;
subplot(1,2,1);
plot(dYearsvec(9:13),DlogTFP_hor_GZsym_DRS(9:13),'-.','LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogTFP_GZsym_DRS(9:13),'LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogTFP_GZ_DRS(9:13),':','LineWidth',3.5,'Color',[.4 .4 .4]); axis tight; grid on; hold on;
ylabel('\DeltalogA');
legend('Horizontal, Symmetric \phi','IO Network, Symmetric \phi','IO Network, Asymmetric \phi','Location','southwest');
title('\DeltalogA with DRS');
subplot(1,2,2);
plot(dYearsvec(9:13),DlogTFP_hor_GZsym_CRS(9:13),'-.','LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogTFP_GZsym_CRS(9:13),'LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogTFP_GZ_CRS(9:13),':','LineWidth',3.5,'Color',[.4 .4 .4]); axis tight; grid on; hold on;
ylabel('\DeltalogA');
title('\DeltalogA with CRS');



figure;
subplot(1,2,1);
plot(dYearsvec(9:13),DlogLWedge_hor_GZsym_DRS(9:13),'-.','LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogLWedge_GZsym_DRS(9:13),'LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogLWedge_GZ_DRS(9:13),':','LineWidth',3.5,'Color',[.4 .4 .4]); axis tight; grid on; hold on;
ylabel('\Deltalog\Lambda');
legend('Horizontal, Symmetric \phi','IO Network, Symmetric \phi','IO Network, Asymmetric \phi','Location','southwest');
title('\Deltalog\Lambda with DRS');
subplot(1,2,2);
plot(dYearsvec(9:13),DlogLWedge_hor_GZsym_CRS(9:13),'-.','LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogLWedge_GZsym_CRS(9:13),'LineWidth',3.5,'Color',[0 0 0]); axis tight; grid on; hold on;
plot(dYearsvec(9:13),DlogLWedge_GZ_CRS(9:13),':','LineWidth',3.5,'Color',[.4 .4 .4]); axis tight; grid on; hold on;
ylabel('\Deltalog\Lambda');
title('\Deltalog\Lambda with CRS');



%% Model-implied STATISTICS
% Table in main text for network multiplier


GZ_decomp2009=[DlogTFP_hor_GZsym_DRS(12)   DlogTFP_GZsym_DRS(12)          DlogTFP_GZ_DRS(12)          NaN                                             ;...
    DlogTFP_hor_GZsym_CRS(12)          DlogTFP_GZsym_CRS(12)      DlogTFP_GZ_CRS(12)      NaN                                             ;...
    DlogLWedge_hor_GZsym_DRS(12)            DlogLWedge_GZsym_DRS(12)        DlogLWedge_GZ_DRS(12)        DlogLWedge_GZsym_DRS(12)/DlogLWedge_hor_GZsym_DRS(12)    ;...
    DlogLWedge_hor_GZsym_CRS(12)        DlogLWedge_GZsym_CRS(12)    DlogLWedge_GZ_CRS(12)    DlogLWedge_GZsym_CRS(12)/DlogLWedge_hor_GZsym_CRS(12)];

disp('Table in the main text');
disp('The Labor Wedge Network Multiplier during the 2008-2009 Financial Crisis');
GZ_decomp2009'





%% Aggregate Labor and Consumption in the data
% aggregate labor in the data
aggLabor    =sum(l_it_mat);
log_aggL    =log(aggLabor);

% aggregate consumption in the data: total consumption divided by aggregate price level
aggPC          =sum(pc_it_mat);
log_idealP     =transpose(diag((v_it_mat')*logp_mat));
log_aggC       =log(aggPC)-log_idealP;

Dlog_aggC     =diff(log_aggC);
Dlog_aggL     =diff(log_aggL);


%% Calculate EQUILIBRIUM OUTPUT AND LABOR
% Parameter Values
epsilon=.5;
gamma=.1;

% Aggregate Labor Wedge in the Data
Dlog_dataLW  =(1+epsilon)*Dlog_aggL-(1-gamma)*Dlog_aggC;

denom_DRS       =(1+epsilon)-etabarDRS*(1-gamma);
denom_CRS       =(1+epsilon)-etabarCRS*(1-gamma);

Gamma_l_DRS     =(1-gamma)/denom_DRS;
Gamma_c_DRS     =(1+epsilon)/denom_DRS;
Lambda_l_DRS      =1/denom_DRS;
Lambda_c_DRS      =etabarDRS/denom_DRS;

Gamma_l_CRS     =(1-gamma)/denom_CRS;
Gamma_c_CRS     =(1+epsilon)/denom_CRS;
Lambda_l_CRS      =1/denom_CRS;
Lambda_c_CRS      =etabarCRS/denom_CRS;

% Equil Labor
% Baseline with GZ
DlogEQLGZ               =Gamma_l_DRS*DlogTFP_GZ_DRS +Lambda_l_DRS*DlogLWedge_GZ_DRS;
DlogEQL_GZsym           =Gamma_l_DRS*DlogTFP_GZsym_DRS +Lambda_l_DRS*DlogLWedge_GZsym_DRS;
DlogEQL_hor_GZsym       =Gamma_l_DRS*DlogTFP_hor_GZsym_DRS +Lambda_l_DRS*DlogLWedge_hor_GZsym_DRS;
DlogEQL_hor_GZasym      =Gamma_l_DRS*DlogTFP_hor_GZasym_DRS +Lambda_l_DRS*DlogLWedge_hor_GZasym_DRS;
% CRS Baseline with GZ
DlogEQLGZ_CRS           =Gamma_l_CRS*DlogTFP_GZ_CRS +Lambda_l_CRS*DlogLWedge_GZ_CRS;
DlogEQL_GZsym_CRS       =Gamma_l_CRS*DlogTFP_GZsym_CRS +Lambda_l_CRS*DlogLWedge_GZsym_CRS;
DlogEQL_hor_GZsym_CRS   =Gamma_l_CRS*DlogTFP_hor_GZsym_CRS +Lambda_l_CRS*DlogLWedge_hor_GZsym_CRS;
DlogEQL_hor_GZasym_CRS  =Gamma_l_CRS*DlogTFP_hor_GZasym_CRS +Lambda_l_CRS*DlogLWedge_hor_GZasym_CRS;


% Equil Consumption
% Baseline with GZ
DlogEQCGZ               =Gamma_c_DRS*DlogTFP_GZ_DRS +Lambda_c_DRS*DlogLWedge_GZ_DRS;
DlogEQC_GZsym           =Gamma_c_DRS*DlogTFP_GZsym_DRS +Lambda_c_DRS*DlogLWedge_GZsym_DRS;
DlogEQC_hor_GZsym       =Gamma_c_DRS*DlogTFP_hor_GZsym_DRS +Lambda_c_DRS*DlogLWedge_hor_GZsym_DRS;
DlogEQC_hor_GZasym      =Gamma_c_DRS*DlogTFP_hor_GZasym_DRS +Lambda_c_DRS*DlogLWedge_hor_GZasym_DRS;
% CRS Baseline with GZ
DlogEQCGZ_CRS           =Gamma_c_CRS*DlogTFP_GZ_CRS +Lambda_l_CRS*DlogLWedge_GZ_CRS;
DlogEQC_GZsym_CRS       =Gamma_c_CRS*DlogTFP_GZsym_CRS +Lambda_l_CRS*DlogLWedge_GZsym_CRS;
DlogEQC_hor_GZsym_CRS   =Gamma_c_CRS*DlogTFP_hor_GZsym_CRS +Lambda_l_CRS*DlogLWedge_hor_GZsym_CRS;
DlogEQC_hor_GZasym_CRS  =Gamma_c_CRS*DlogTFP_hor_GZasym_CRS +Lambda_l_CRS*DlogLWedge_hor_GZasym_CRS;



%% FIGURES for Model-Implied TFP and Labor Wedge
% Figures in online appendix

figure;
plot(dYearsvec,DlogTFP_GZ_DRS,'-.','LineWidth',3.5,'Color',[25/255 25/255 122/255]); axis tight; grid on; hold on;
plot(dYearsvec,DlogTFP_GZ_CRS,'LineWidth',3.5,'Color',[128/255 0/255 0/255]); axis tight; grid on; hold on;
%plot(dYearsvec,Dlog_dataTFP,':','LineWidth',3.5,'Color',[0 128/255 128/255]); axis tight; grid on; hold on;
ylabel('\DeltalogTFP');
legend('\DeltalogTFP, DRS','\DeltalogTFP, CRS','Location','southwest');  

figure;
plot(dYearsvec,DlogLWedge_GZ_DRS,'-.','LineWidth',3.5,'Color',[25/255 25/255 122/255]); axis tight; grid on; hold on;
plot(dYearsvec,DlogLWedge_GZ_CRS,'LineWidth',3.5,'Color',[128/255 0/255 0/255]); axis tight; grid on; hold on;
plot(dYearsvec,Dlog_dataLW,':','LineWidth',3.5,'Color',[0 128/255 128/255]); axis tight; grid on; hold on;
ylabel('\Deltalog\Lambda');
legend('\Deltalog\Lambda, DRS','\Deltalog\Lambda, CRS','\Deltalog\Lambda, data','Location','southwest');



%% Decomposition of Output and Labor by TFP and Labor Wedge
% Figure in online appendix

figure;

subplot(2,2,1);
Y=[Gamma_c_DRS*DlogTFP_GZ_DRS' Lambda_c_DRS*DlogLWedge_GZ_DRS'];
h = area(dYearsvec,Y); axis tight; grid on; hold on; 
h(1).FaceColor = [148/255 0/255 211/255];
h(2).FaceColor = [238/255 130/255 238/255];
legend('TFP','Labor Wedge','Location','southwest');  
ylabel('\DeltalogC');
title('\DeltalogC, DRS');

subplot(2,2,2);
Y=[Gamma_c_CRS*DlogTFP_GZ_DRS' Lambda_c_CRS*DlogLWedge_GZ_DRS'];
h = area(dYearsvec,Y); axis tight; grid on; hold on; 
h(1).FaceColor = [148/255 0/255 211/255];
h(2).FaceColor = [238/255 130/255 238/255];
ylabel('\DeltalogC');
title('\DeltalogC, CRS');

subplot(2,2,3);
Y=[Gamma_l_DRS*DlogTFP_GZ_DRS' Lambda_l_DRS*DlogLWedge_GZ_DRS'];
h = area(dYearsvec,Y); axis tight; grid on; hold on; 
h(1).FaceColor = [148/255 0/255 211/255];
h(2).FaceColor = [238/255 130/255 238/255];
ylabel('\DeltalogL');
title('\DeltalogL, DRS');

subplot(2,2,4);
Y=[Gamma_l_CRS*DlogTFP_GZ_DRS' Lambda_l_CRS*DlogLWedge_GZ_DRS'];
h = area(dYearsvec,Y); axis tight; grid on; hold on; 
h(1).FaceColor = [148/255 0/255 211/255];
h(2).FaceColor = [238/255 130/255 238/255];
ylabel('\DeltalogL');
title('\DeltalogL, CRS');

disp('decomposition percentage: TFP vs Labor Wedge')
PctCons2009_GZ_DRS=(Lambda_c_DRS*DlogLWedge_GZ_DRS(12))/(Gamma_c_DRS*DlogTFP_GZ_DRS(12)+Lambda_c_DRS*DlogLWedge_GZ_DRS(12))
Pctlabor2009_GZ_DRS=(Lambda_l_DRS*DlogLWedge_GZ_DRS(12))/(Gamma_l_DRS*DlogTFP_GZ_DRS(12)+Lambda_l_DRS*DlogLWedge_GZ_DRS(12))

PctCons2009_GZ_CRS=(Lambda_c_CRS*DlogLWedge_GZ_DRS(12))/(Gamma_c_CRS*DlogTFP_GZ_DRS(12)+Lambda_c_CRS*DlogLWedge_GZ_DRS(12))
Pctlabor2009_GZ_CRS=(Lambda_l_CRS*DlogLWedge_GZ_DRS(12))/(Gamma_l_CRS*DlogTFP_GZ_DRS(12)+Lambda_l_CRS*DlogLWedge_GZ_DRS(12))

